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Highly frustrated spin systems represent a central and challenging problem in condensed mater 
physics. To this problem, we introduce an algorithm based on mixed projected entangled pair states 
(ra-PEPS), which is a novel type of tensor network. We use the famous Kitaev model on an infinite 
honeycomb lattice, which can be solved exactly, as a benchmark. With very limited parameters and 
finite scaling, our calculation results are in good agreement with the exact results, indicating the 
efficiency of our algorithm. After presenting the benchmark, we investigate the Kitaev-Heisenberg 
model, which was proposed to describe the effective magnetic momentum interaction in iridate 
Na2ir03 which may be used to realize the spin liquid phase. However, our calculations suggest that 
the gap less spin liquid phase is not robust at the thermodynamic limit, and thus this phenomenon 
is very difficult to observe. 
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Frustrated spin systems play an essential role in two- 
dimensional many-body physics and have attracted the 
attention of condensed matter physicists for decades [2]. 
In particular, some exotic ground states, such as the 
quantum spin liquid state (which has recently been 
strongly supported experimentally in the Kagome lattice 
by Y. Lee [3] et a/), may play a key role in explaining high- 
T c superconductivity, as proposed by P. W. Anderson[l]. 
The Kitaev model on a honeycomb lattice [4], which can 
be solved exactly and supports a spin liquid phase, is an 
important frustrated model in two-dimensional physics. 
In addition, this proposed model was proposed that it 
may be realized in a cold atom system [5], or in the 
iridate Na2lr03 [19]. The exact solution shows that the 
spin liquid phase is gapless, and a gap opens when a 
magnetic field, which breaks the time-reversal symmetry, 
is applied to this model as a perturbation. Unfortunately, 
general two-dimensional frustrated spin models - even 
slightly modified Kitaev models such as the Kitaev model 
perturbed with a magnetic field - cannot be solved 
exactly. 

Due to the lack of analytical methods for general 
many-body systems, numerical analysis remains 
the main method for understanding such physics. 
However, traditional algorithms, such as the quantum 
Monte Carlo (QMC) method and density matrix 
renormalization group (DMRG) method, have failed 
to simulate frustrated systems in two dimensions: 
QMC suffers from the notorious 'sign problem' [15] 
while DMRG is limited to one dimension and generally 
lose its power in higher dimensional systems [6, 7]. 
Fortunately, the recently developed tensor network 
(TN) algorithms, such as, the algorithm based on 
projected entangled pair states (PEPS) [16-18] which is a 
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natural generalization of the DMRG algorithm to higher 
dimensions, has shown great potential for addressing 
this problem. The algorithm based on infinite PEPS 
has been used to obtain a reasonable phase diagram 
of the frustrated antiferromagnetic J\ — J2 Heisenberg 
model on a checkerboard lattice[14]. More recently, 
a spin liquid phase was claimed near the maximally 
frustrated region (J2 ~ 0.5Ji) for the spin 1/2 Ji — J2 
antiferromagnetic Heisenberg model on a square lattice 
for the TN algorithm[12, 13]. These results show that 
the TN algorithms are powerful tools for exploring the 
frustrated systems and for evaluating some new physics 
beyond the traditional methods. 

Despite the success of the TN algorithm, the study of 
highly frustrated large systems, such as those represented 
by the Kitaev model and related models, remains a 
challenge. The general TN algorithm, such as the 
algorithm based on PEPS, cannot give the appropriate 
results for the Kitaev model, because the limited ability 
of the calculation constrains the bond dimension D to 
a relatively small number (typically < 10). Thus new 
methods must be developed to highly frustrated systems 
with the current calculation abilities. In this letter, we 
introduce an algorithm based on infinite mPEPS, which 
is a mixture of the bosonic PEPS (bPEPS)[17] and the 
fermionic PEPS (fPEPS)[8, 9]. We will demostrate that 
this novel method can be used to efficiently study highly 
frustrated systems with a relatively small parameter D. 
We use the Kitaev model as a benchmark and obtain 
satisfactory results, for both the energy, and the corre- 
lations of this model. After the benchmark calculation, 
we apply this method to the Kitaev-Heisenberg model, 
which is proposed to describe the effective magnetic 
momentum interaction model of the iridate Na2lr03 
[10, 19], to obtain a phase diagram. Our calculations 
indicate that the gapless spin liquid region is very narrow 
in the thermodynamic limit and that would thus be very 
difficult to observe in the iridate Na2lr03, even with the 
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Kitaev -Heisenberg model. 

Mixed Projected Entangled-Pair States The TN 

algorithm is a variational method based on TN states 
with special structures. A general TN state is defined 
in two steps: first, we describe a configuration of 
entangled states between virtual particles (we call each 
entangled state a bond, connected by a line, as shown 
in Fig.l), which determines the structure and the upper 
bound of the entanglement of a TN state; second, we 
define projectors to project the virtual-particle space into 
real physical space. The different TN algorithms are 
distinguished by the structure of the TN state, such 
as PEPS, or the string bond state. The variational 
parameter space of TN algorithms is determined by 
the parameters of the projectors, the number of which 
is dependent on a polynomial of the dimension of the 
bond D, the dimension of the physical space d and the 
number of projectors (as determined by the structure 
of the lattice and the symmetry of the state). For the 
special TN state, PEPS, the configuration of the bonds 
is completely determined by the lattice where each edge 
of the lattice corresponds to a bond, and each site has 
an on-site projector. Without a loss of generality, we use 
the honeycomb lattice as an example. In general, the 
entangled pair can be written as: EPR e = |z) ei |i) e2 , 
where \i) is the state of virtual particle with dimension 
D and e denotes an edge connecting sites e\ and e2 in 
the lattice. Algebraically, the state \i) can be represented 
by local creation operators from the vacuum through the 
second quantization method. In the original PEPS case, 
the entangled pair state can be represented as bosonic 
operator a) as: EPR e = ^ a t\ a e\ \V ac )- 

Similarly, the projector of each site on a honeycomb 
lattice can also be represented as bosonic operators 
as: plhh] = £ pi) . k T^b^aia{a k z \vac){Vac\, where 
[hfa] is the coordinate of the site, is the bosonic 
creation operator of the physical particle, and a Xl a yi a z 
are the bosonic annihilation operators of the virtual 
particles connected to edges x,y and z, respectively 
(see Fig.l). In this case \vac) is the vacuum of the 
physical space, and \Vac) is the vacuum of the virtual 
space. Using the representation of the entangled pairs 
and projectors, the PEPS state can be represented as: 

= hi u i 2 p[hh] rL EPR e> where a11 of the virtual 
particles are projected as physical particles. 

In general, the virtual particles in an EPR pair 
are not necessarily bosons. In solving the many-body 
physics of fermionic systems, it is natural to extend 
the bPEPS formulism to the fPEPS formulism. This 
can be completed by replacing the bosonic operators in 
the entangled pairs and projectors, as shown in Eq.(l) 
with fermionic operators. In the fPEPS formulism, 
the parity of the projectors at different sites of the 
lattice, which determine the exchange character of the 
projectors, should be the same even. This requirement 



guarantees that the fPEPS is well defined, that is, the 
order of the projectors is not important for the state 
beyond the global phase. 

Here we focus on a frustrated spin system, which 
exhibits bosonic statistics in physical space. Thus, 
we propose an algorithm based on the mixed PEPS 
(mPEPS), which introduces fermionic statistics to the 
virtual-particle space (same as fPEPS) and bosonic 
statistics in the physical-particle space (different from 
fPEPS). Therefore, the projectors of the mPEPS 
should be expressed in the second quantized form as: 
^ P ,i,j,k T p ij k b^ p ala J y a^\vac) (Vac\, where 6 1 " satisfies: 
= Sij and a x ,a y ,a z satisfy: {a*, at} = Sij. To 
validate this definition, the parity of the fermionic 
operators in the projector should also be same for every 
site (for example, it could be even). Although there is 
some overlap, the state spaces of the PEPS, fPEPS and 
the mPEPS are expected to be different, even for the 
same parameters. 

After the defining mPEPS, the ground state of a sys- 
tem in variational space can be found by the imaginary 
evolution method which is performed in the same manner 
as in the PEPS method; in particularly, we use the simple 
update method to approximate the environment in the 
same manner as PEPS in two dimensions [24]. When 
the ground state is stable under the imaginary evolution, 
the physical quantities of the system can be calculated 
through a complex contraction process (for more details, 
see the Supplementary Information). In the contraction, 
the bosonic and fermionic operators must be carefully 
handled(see in the Supplementary Information). 




FIG. 1: (Color on line) Representation of the mPEPS on 
a honeycomb lattice: each line denotes an entangled pair 
state; the blue balls denote virtual fermions; and the red balls 
denote bosons; the circle around each site represents an on- 
site projector P. Different link directions are denoted by x,y 
and z, which are used in the Kitaev model. In the following 
calculations, we use an infinite lattice with translational 
symmetry and six different projectors, Pi, P2, • • • , Pq 
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FIG. 2: (Color on line) Nearest Neighbor (NN) correlation 
(a z (j z )z of the Kitaev model along the line J x = J y , which 
is indicated in the left upper inset by a dashed line. The red 
dots represent the results obtained by the mPEPS method 
with finite scaling of D; the black line shows the exact 
results of (a z a z ) which is given as the formula (afa z ) = 

J BZ cosO (fci, k2)dk\dk2 [11]; the blue triangles display 
the results obtained by the PEPS method using the same 
parameters of mPEPS. The right lower inset shows the (§^=r}, 
which displays the transition point near 0.46. 



Kitaev Model as the Benchmark Here, we use 
the Kitaev model, which is a highly frustrated spin 
model in two dimensions, as a benchmark of the mPEPS 
algorithm. The Kitaev model is defined on an in- 
finite honeycomb lattice [4], and its interaction along 
different directions (see Fig.l) vary strongly as: H = 
E 7 -/mfc J-y a i°'] - This m °del can be exactly solved with 
two resultant phases: the gapped phase and the gapless 
spin liquid. There are two types of quasi-particles in 
the gapped phase: fermions and vortices. The vortices 
associated with a Z2 gauge field and a factor —1 will 
appear when the fermion moves around the vortex. A gap 
opens when a magnetic field is applied to the gapless spin 
liquid phase. The non-abelian anyon, which may play 
an important role in topological quantum computation, 
will also appear in this phase. In our calculation, we 
normally set the parameters of the model to satisfy J x + 
Jy + Jz = 1 and J x — J y , which gives a line connecting 
a vertex to the middle point of the corresponding edge 
(see the red dashed line in the inset of Fig. 2). After 
conducting the stability tests of the algorithm with 
respect to the parameters (particularly, the truncation 
parameter D cut in the contraction process, see details 
in the Supplementary Information), we compare the 
energy of the system along the chosen line with the 
exact solution. Result show good agreement (see the 
Supplementary Information). In addition, the nearest 
neighbor (NN) correlation (<j z <j z ) which corresponds to 



(§Jz) a l° n g the line is also shown in the Fig. 2. The black 
line shows the exact value of this correlation as derived 
in[ll], which is in good agreement with our calculation 
(red dots). The inset shows the (Jj^), which exhibits a 
sharp transition point near 0.46 and indicates a phase 
transition in this model, although it slightly deviates 
from the exact result of 0.5. 

We must emphasize that the virtual dimension of the 
system D is rather limited by our calculation ability. 
We can only calculate results for the dimension up to 
8. And we conduct the finite scaling for this limited 
dimension in this model, and the points shown in Fig. 2 
were obtained with those finite scaling. To compare our 
findings with the traditional PEPS case, we also present 
the results from the PEPS with the same parameters 
and finite scaling employed in Fig. 2 (blue triangles). It is 
clear that the present method improves the results obtain 
with limited calculation ability. More importantly, the 
PEPS method does not give an indication of the phase 
transition, while our mPEPS technique can. However, in 
simpler cases, such as two-dimensional Ising model with a 
transverse magnetic field on the honeycomb lattice, the 
mPEPS and PEPS methods give nearly the equivalent 
results. 

Analytical result [4] suggest that the gapless spin liquid 
will exhibit a gap upon the addition of a magnetic field 
in the direction e x + e y + e z ((1,1,1) direction) as a 
perturbation, which breaks the time-reversal symmetry. 
Unfortunately, the Kitaev model with magnetic field can- 
not be solved exactly. We can only study this transition 
by numerical methods. Our mPEPS calculations show 
this transition directly by — where h is the magnitude 
of the magnetic field (see Fig. 3), and these calculations 
strongly suggest a first order phase transition occurring 
at very low magnetic fields, h < 0.005. 

Application to the Kitaev-Heisenberg Model 

Having verified the mPEPS method in the Kitaev model, 
we extend the application of this method to the Kitaev- 
Heisenberg model. The Hamiltonian of this model can 
be written as: Hhk — (1 — &)Hh + IolHk, where Hh is 
the standard Heisenberg interaction: V\ „, „ (J 1 a 1 and 
the Hk = ^l* 7 ] i s Kitaev interaction. With the 

Khaliullin transformation [19], which performs rotations 
on different sites, this model can be exactly solved at 
the point a — 0.5 and a = 1. When a = 0.5, the 
ground state of the transformed system can be exactly 
solved, and is found to be a ferromagnetic state. With 
the inverse Khaliullin transformation, the ferromagnetic 
state will be transformed to the stripy anti-ferromagnetic 
state. When a = 1, this model is equivalent to the Kitaev 
model, with J x = J y = J Z1 whose ground state can be 
exactly solved as a gapless spin liquid state. For the 
parameters beyond a = 0.5 and a — 1, the model can not 
be exactly solved and can only be examined by numerical 
calculations. Former studies have employed several 
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numerical methods, such as exact diagonalization[19, 21] 
and DMRG[20], and have found three different phases: 
the Neel order phase, the stripy anti-ferromagnetic phase 
and the spin liquid phase. These previous calculations, 
which were based on the systems of relatively small finite 
size[19] or quasi one-dimensional [20], showed that phase 
transition points occur near a = 0.4 and a = 0.8 at zero 
temperature. 

We applied the mPEPS method to this model on an 
infinite two-dimensional lattice which naturally satisfies 
the thermodynamic limit. We directly calculate the spin 
configuration of the ground state with a = 0.5 which is 
indeed a ferromagnetic state, as predicted by the exact 
results. As with the benchmark calculation of the Kitaev 
model with J x = J y = J z = 1/3 (corresponding to a = 
1), here we also give which can be used to determine 
the phase transition point, as shown in Fig. 4. 
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FIG. 3: (Color on line) The transition according to the 
modulation of — ^ with a magnetic field h. The Kitaev 
model is initially in the spin liquid phase with parameters 
J x = Jy = Jz = 1/3, as denoted by the green point. 
The red dots indicate results for the parameter D — 4, and 
the black squares display the results for D — 6. Although 
there is numerical difference between the two parameters 
(without finite scaling for D), the transition is sufficiently 
sharp enough and at the same position. This finding provides 
strong evidence that there is a first phase transition occurs 
near h = 0.005 under finite scaling. 

In our infinite system with six-site periodicity, the 
second transition point is dramatically different from 
the former result (a = 0.8), while the first transition 
point, at a = 0.46, is slightly different from a = 0.4. 
To confirm our calculation, we also determine the NN 
correlation (<?<?), which is used as an indication of the 
phase transition in reference [19]. The results also support 
the occurrence of a phase transition at the same position. 
In our calculations, the second phase transition point 
occurs at a « 0.98 for D = 4, and at a « 0.99 
for D = 6. This findings suggest that the transition 
point will tend to a = 1 with the finite scaling for D. 
Our results indicate that the gapless spin liquid phase 



will be destroyed by the Heisenberg perturbation very 
quickly, and thus this phase is not robust for this type 
(Heisenberg) of perturbation. The direct implication 
of our results indicates that it would be very difficult 
to observe a Kitaev-type gapless spin liquid in iridate 
Na2lr03. The findings of a decreasing spin liquid phase 
in the infinite system (at the thermodynamics limit), 
compared to the finite exact diagonalization and DMRG 
in quasi one-dimension with cylinder boundary condition, 
is not so surprising. The size of the system plays a 
subtle role in the spin liquid phase. Z. Meng et al have 
reported [22] a spin liquid phase between the semimetal 
phase and an antiferromagnetic Mott insulator in the 
Fermi-Hubbard model on the honeycomb lattice, with 
the inclusion of 648 sites, by using projective QMC 
method. However, when S. Sorella et al [23] reconsidered 
this problem by including up to 2592 sites, the evidence 
of spin liquid was lost. In our case, the infinite lattice 
may play a similar role. 
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FIG. 4: (Color online) The values of g (a) and the NN 
correlation (aa) (b) versus a for D = 4 (black and green 
squares) and D = 6 (red and blue dots). The left lower inset 
shows the magnification near a = 1 and the right upper inset 
displays the configuration of the spin under the Khaliullin 
transformation for a = 0.5. It is clear that two phase 
transitions occur. The first transition point is near a — 0.46 
which is only slightly different with 0.4. However, the second 
phase transition occurs at a ~ 0.98 which is qualitatively 
different from 0.8. Our results show that the spin liquid region 
is very narrow and that the gapless spin liquid will be rapidly 
destroyed by the Heisenberg perturbation. 

To summarize, we have introduced a novel type of 
TN (mPEPS) algorithm to attack the highly frustrated 
spin systems. Our method is successful in simulating the 
Kitaev model, with respect to both the energy and the 
correlations of the ground state. Having established this 
method, we applied it to the Kitaev-Heisenberg model. 
Our calculations show that the spin liquid region is very 



narrow, and thus this state would be very difficult to 
observe experimentally. The success of our method opens 
the possibility of considering virtual particles in a tensor 
network with more complex statistical characteristics, 
such as the fractional statistics of anyons. These new 
statistics may be convenient for studying the system 
which excitations have the similar statistics. 
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Nos. WK2470000006, WK2470000004, WJ2470000007 
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Supplementary 
Imaginary Time Evolution 

After chosen the variational space of the state, i.e., 
mPEPS with given parameter D and d, we use the 
imaginary time evolution method[16-18] to find the 
ground state. The imaginary time evolution could be 
summarized as the following formula: 



The ) is the initial state, which can be chosen 
randomly, after long enough imaginary time evolution, 
the final normalized state will be close to the ground 
state with well controlled precision, as long as the initial 
state is not orthogonal with the ground state (with zero 
measurement to chosen these states). Generally, we can 
only deal with the local evolution operators (the time 
evolution of a system is a global operator), and the 
Trotter expansion can replace the global operators by 
local operators with controllable precision[16-18]. For 
example, in Kitaev model, Hamiltonian H could be split 
as H = H x + H y + H z , where Hi,(i = x,y,z) is the 
part of Hamiltonian acting on i-bond. With the help of 
Trotter expansion, we could split the evolution operator 
in second order as: 

e -Hdt _ e -(H x +H y +H z )dt 

_ e -H x dt e -H y dt e -2H z dt e -H y dt e -H x dt _|_ 0(d£ 3 ) 




(d) (c) 



FIG. Al: SVD in imaginary time evolution: (a) Act the 
evolution operator e~ dtH (the parallelogram blue tensor) on 
two contracted tensors. Then we obtain a bigger tensor as 
indicated by (b). (c) Use SVD to cut the big tensor and (d) 
reserve only D biggest singular values to make the tensors 
back to their original dimension. For convenience, before any 
contract processes of the tensors, we contract the diagonized 
matrix, which is the approximation of the environment, into 
the nearest tensors. 

Contract 

After the imaginary time evolution, we suppose that 
we have found the ground state. Then the key task is 
to find the values of some physical quantities which can 
be compared with the experiment. This task is realized 
by the contract processes, which including the contract 
of tensors and entangled pairs. The physical quantities 
in a given ground state |^) can be expressed as 



With the Trotter expansion, it seems that we just need 
to act the local evolution operators on the initial state 
with the given order and obtain the ground state finally. 
However, the dimension of the bond D will increase 
exponentially with the evolution time t, if local operators 
are applied without approximation. Therefore, as in the 
general TN algorithm, approximations should be made 
after each step of evolution. One of the simplest way 
to cut the tensors back to the initial parameter space is 
the singular value decomposition (SVD), which maintain 
the important part of the tensor (see Fig.Al). It is 
straightforward to use SVD in the case of PEPS, while 
it needs extra operations in the case of mPEPS (also 
in fPEPS) with fermionic statistics. Since the tensors 
in mPEPS are set with even parity, the matrix in the 
SVD process can be re-indexed as block diagonal matrix 
according to parity of the tensor's indices (one block with 
odd parity and the other with even parity), then we can 
apply the SVD and cut off in each of the block. 



For simplicity, we use (^ g \^ g } as an example to show 
how to contract to calculate the physical quantities. The 
calculation of physical quantities is similar. Suppose we 
have a 2D honeycomb lattice and each site is labeled with 
two indices i,j. We have a state in the form of projected 
entangled-pair State: 

w=u pM n epi wm' 

ij iijii2j2 

where EPR [ilil][i2j2] = J2 V \ v )[iih] \ v )[iih] and [hji]^} 
denotes the edge e, linking site [iiji] and site [22.72] ■ 
Normally, translational invariance of tensors will be 
imposed in infinite lattice. Actually, we have the 
corresponding bra vector as 

m= n EpR U]^]ii^ ]t ' 
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where EPRf.^.^.^ = E^k^Hw The Product 
order is reversed which will make a difference in the 
contract procedure below. The scalar product of bra 
vector and ket vector is: 

W)= II EPR Um^U pm 

lljllljl ij 

u p[ij] n Ep^iiii^] 

ij UJl*2j2 

Because we are considering mPEPS where tensors have 
even parity, it is free to change the orders. We would like 
to rearrange the order as: 

im= n EPR U^] 

iijiiiji 

Y[ P m P M Yl EPR [ilil][22i2] 

ij i\j\iiji 

The product of pMlpM w m contract over physical 
indices and we symbolize it as K^. 

Kij P^pM 

= Yl T pvllyV z \ V Z V y V x) (P\pV X VyV z \P) (V X VyV z \ 
= E T m^ T PvlvyV z I W^} (V X V V V Z | 

Y K ^ ] v X Vy-v y v^v z l^} K) ( V V I , 

where repeated indices mean contraction. Notice that 
K Wo x Vy:v y v- z v z is not J ust the contraction of T%jj*^ 
and TpJ] v Vz , but with extra minus signs resulting from 
commutation of operators. 

Therefore, the scalar product is 

<#/,)= n EPRt i . i][i2 . 2] n^- n EPR [ilMi2j2] 

i\j\iiji ij i\j\iiji 

What we should do next is to contract over all virtual 
indices, ie. to product all the K tensors by EPRs. 
Even though we consider here is a infinite lattice, we 
suppose there are boundaries in very remote areas. We 
contract row by row many times from boundaries until 
the boundary tensors are stable. After we contract the 
boundary with a row of tensors, the boundary tensors 
will get thicker, like the case of imaginary time evolution, 
and we need to do approximations to cut them back to 
their initial shape. Normally we use a specific form of 
boundary tensors to make the approximation, for the 
sake of simplicity (see Fig.A2 for details). 




FIG. A2: contract procedures: (1) Contract boundary 
tensors with one row of K tensors; (2) trace the new row 
of tensors with their conjugate; (3) Seek left and right 
environment tensors by iterative product; (4) Obtain L and R 
tensor by cutting left and right environment using eigenvalue 
decomposition, Contract L, the fixed tensor, and R together; 
(5) Use SVD to make the approximation; (6) Contract L's 
inverse and R's inverse tensors respectively. 

The main tool used to make this approximation is 
still SVD: We fix two tensors, and trace over left and 
right parts of boundary tensors to gain the left and 
right environments of the two fixed tensors. Now we 
contract over left environment tensor, fixed tensors, right 
environment tensors together to become a bigger tensor. 
We use SVD to cut the bigger tensor, then cancel the 
environment effect by contracting inverse tensors of the 
left and right environments. It is extremely obvious if you 
write down the approximation judgement explicitly: 

min(Ti(B - Bo)HB - Bo)) 

= min Tt((G l U[U^G r - G L GG R )\G L U[U' 2 G R - G L GG R )) 
u[u' 2 

= min (Tr (G R U? U? G{G L U[U^G R ) 
u 1 u 2 

-Ti(G i R U?U' 1 t G t L GLGG R )Tr(G f R G t GlGLUiu!>G R ) + 
Ty(G r G ] GIG l GG r )) 

= min (Tr(i? t t r 2 t t r ( t if/^2)Tr(7? t f/2 t t r { t LG) 

u x u 2 

-Ti{R ] G ] LU[U' 2 ) + Tr(R*G*LG)) , 

where Bq is the old boundary contracting with a row 
of tensors, B is the new boundary with only two tensors 
different from the old boundary. G is the unit cell of Bq. 
Gl Gr are perspectively the left and right part of G in 
the row of £?o. L is just G^ l Gl and R is just G^ r Gr. 
The above equation is just binary form of U[ and U^. 
Obviously, Tr(£ - B )i(B - B ) reaches its minimum 
when best approximates \^RJgVL. See 

also Fig.A3. 

This example of (^1^) showed our method for contract. 
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It is similar to calculate average value of local operators 
(ip\0\ip), except that the sites where operators lie should 
be kept until final contract. 

There is one more point worth mentioning: When 
we use mPEPS, we choose the virtual annihilation and 
creation operators to be anti-commutative. So we must 
be extremely careful about the minus sign when we 
contract over virtual indices. Similarly and more easily 
to overlook, in SVD, minus signs also appear when we 
change the order of virtual indices in order to obtain 
the original shape of tensors. It could be explained as 
equations below. We take a simplified version as an 
example. T is a big tensor which needs to be split by 
SVD. 



U x U 2 U x U 2 Ui u 2 




f 



U'j U' 2 U'j U' 2 U'j u' 2 



V 



Ui U 2 U'i U' 2 Ui u 2 



B 



T- - 

-J- uu,vv 

ETJ- 9- V' 
^ uu,ss^ ss v SS,VV 

ss 

— ^ ^ Uuu,ss Vss^VV 

ss 

— ^m,sTsl^sl,s2s2^2s2,i;i; 
slsls2s2 

= EPR 1 E^ >iTsl |i2il) (al^la.a^JEPR 

sTsls2s2 

= EPR f (_^ ^ u ,,T sl (-) (p(iT)+p(sl))p(i5) 

slsls2s2 

|5I)<al||i2>(a2|^3. 2i0t ,)EPR 

= EPR f (_^ f/^,jT sl (-) (p(5T)+p(al))p(5T) 

slsls2s2 

|iT)( S l||52>( S 2|y i5s2j ,jEPR 
= EPR f ( ^ cC,^.i 

slsls2s2 

|iI>( S l||52>< S 2|V^ 2>ii jEPR 

= contract U" tensor with V tensor 



The first equation is normal SVD; the second is 
to absorb singular values; and rest equations are to 
extract EPR from tensors. EPR is ^ s2 = sl \s2)\sl) , 
correspondingly EPR^ = J] s2=sl ( ( sl|( ( s2|, and so as their 
Hermit ian conjugate. In the equations, you can see that 
a minus sign show up. Note that 

p(s2) in the fifth equation is replaced by p(sl) since it 
is not zero only when si = s2 and si = s2. p(*) is 
the function that maps quantum numbers to parity. U 
and V is the tensors that U' and V absorb square of 
the singular values respectively. U" is the tensor after U 
absorbs (-)(p( s_1 )+^i)M^). 



FIG. A3: The scheme of approximation in the language of 
variational method. Use the second row of boundary tensors 
to approximate the first row. Substitute the first row with 
the third row. And keep on this procedures until stable. 



Additional Calculation of the Kitaev Model 

In order to shows the stability of our results, we 
increase the dimension of the truncation n the contract 
from D cut = 3 2 to D cut = 7 2 . We can see in the 
following figure that when D cut > 4 2 with D = 8, the 
NN correlation of the system is relatively stable. 



15 20 25 30 35 40 45 50 55 

Dcut 



FIG. A4: The diagram shows that when D = 8, the NN 
correlations (o- z a z ) are stable as D cut increases from D cut = 
4 2 to D cut = 7 2 . 



With the stable algorithm, we calculate the energy 
of the Kitaev model which is well agree with the exact 
energy 
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— i • 1 • 1 1 1 1 1 1 1— 

0.0 0.2 0.4 0.6 0.8 1.0 

Jz 



FIG. A5: The diagram shows the ground state energy of 
Kitaev model by fPEPS.(ST stands for exact results. FS 
stands for the results after finite scaling of D.) 



Since the limit of our calculation, we should do some 
finite scaling. However, the finite scaling of D is not 
clear and well defined. Here we use linear of polynomial 
fit with l/D, as showed below. 



0.71 5 -, 




0.670 -I 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 

0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 

1/D 

FIG. A6: Finite D Scaling of (a z a z ) z at Jz=0.4. We choose 
it to linear- fit with 1/D. 




FIG. A7: Khaliullin Transformation: The spins on square 
sites are fixed. The spins on circle sites are rotated by tv 
along x bond. The spins on triangle sites are rotated by tt 
along y bond. The spins on pentagon sites are rotated by tv 
along z bond. 



Khaliullin Transformation for Kitaev- Heisenberg 
Model 

This transformation first proposed by Khaliullin and 
Chaloupka not only makes the anti-ferromagnetic state 
to ferromagnetic state, but also makes Kitaev-Heisenberg 
model exactly solvable at a = 1/2. This transformation 
requires different spins in the lattice to rotate about 
different axis. To be specific, we first choose a set of spins 
which are positioned on third nearest neighbor sites at 
opposite corners of the hexagons throughout the lattice, 
and hold these spins fixed. We next rotate the three 
nearest spins by tt about the spin axis corresponding to 
the bond which connects it to the fixed spin. The effect 
of this operation therefore transforms the Heisenberg 
Hamiltonian to a mixture of Heisenberg Hamiltonian and 
Kitaev Model, and Kitaev Hamiltonian invariant. 

Hh — > —H H + 2Hk 
Hk — > Hk 

Thus the Kitaev-Heisenberg model is transformed as: 

H->-(l- a)H H - 4(a - 1/2)H K (Al) 

It is obvious to find out in the transformed Hamilto- 
nian, when a = 1/2, it is just H = — 1/2Hh Heisenberg 
Hamiltonian with a negative coupling coefficient. So 
the state at a = 1/2 after the transformation is a 
ferromagnetic state. We can operate an inverse trans- 
formation to send the ferromagnetic state back to an 
anti- ferromagnetic state. As we have explained before, 
the anti-ferromagnetic state is harder to calculate, since 
stable results are not easily required. So in the letter, 
we just use the Hamiltonian after the transformation to 
discuss the quantum phase diagram of Kitaev-Heisenberg 
model. 

In order to weaken the influence by Heisenberg inter- 
action and make the spin liquid region broader and easier 
to see, we re-parameterize the Kitaev-Heisenberg model 
to: 

H -+ ~ l -^-H H - 4(a' - l/2)H K (A2) 

The new parameter a' and the old one a can be related 
by a simple formula, if Hk is proportional to Hh to 
the same extent. Besides, the new Hamiltonian Eq.A2 
preserves the same behaviors as Eq.Al at a' = 1 and 
a' = 1/2. In the paper, we discuss Eq.A2, and note 
that the spin liquid region of Eq.Al is actually much 
much narrower. The a' region of spin liquids is 0.98 < 
a' < 1, as the part below suggests. Correspondingly, 
the real parameter region of spin liquids before re- 
parameterization are around 0.997(a(l. In the paper, 



we just use the new Hamiltonian Eq.A2 to further our 
discussion. 

The configuration of the ground state of the Kitaev- 
Heisenberg model with a — 0.5 can be calculated by our 
PEPS method as Fig.A8 



FIG. A8: The above figure shows the spin orientation of a = 
1/2 Kitaev-Heisenberg model after Kahliullin transformation. 
And with the inverse transformation, the spin pattern returns 
to stripy anti-ferromagnetic order. 



